■ UN,

Marine Physical Laboratory

APPLICATIONS OF BISPECTRAL ANALYSIS

A. M. Richardson and W. S. Hodgkiss

MPL Technical Memorandum 433

MPL-U-43/93 August 1993

DT1C

electe MUR 1 4 1995

G

Approved for public release.

D-TIC QU'ALHT INSPECTED i

University of California, San Diego Scripps Institution of Oceanography

19950310 041

REPORT DOCUMENTATION PAGE

Form Approved OMB No, 0704-0188

Public reporing burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering ana maintaining the data need ed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204. Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

1. Agency Use Only (Leave Blank).

4. Title and Subtitle.

2. Report Date. August 1993

3. Report Type and Dates Covered. Technical Memorandum

Applications of Bispectral Analysis

5. Funding Numbers. N0001 4-92-K-2005

6. Author(s).

A. M. Richardson and W. S. Hodgkiss

7. Performing Monitoring Agency Names(s) and Address(es).

University of California, San Diego Marine Physical Laboratory Scripps Institution of Oceanography San Diego, California 92152

9. Sponsoring/Monitoring Agency Name(s) and Address(es).

Naval Research Laboratory 4555 Overlook Avenue, S.W.

Washington, D.C. 20375-5000 Code 5160

Project No. Task No.

10. Soonsoring/Monitoring Agency Report Number. 3

12a. Distribution/Availability Statement.

Approved for public release; distribution is unlimited.

12b. Distribution Code.

13. Abstract (Maximum 200 words).

Estimates of the power spectral density have been found to be very useful in a variety of signal processing applications over the last several decades. Higher order spectral contain information not present in the power spectrum and recently, estimates of higher order spectral have been shown to be useful in certain signal processing problems. In particular, estimates of the bispectrum and bicoherence have been found useful in detecting non-Gaussianity and non-linearity, in system identification, and in detecting transient signals.

This paper is concerned with the estimation of the bispectrum and bicoherence of underwater acoustic signals. We first in¬ troduce and describe some of the properties of the bispectrum and bicoherence. The bispectral estimation algorithm as imple¬ mented here at the Marine Physical Laboratory is then described in detail. Finally, several bispectra are estimated from actual data and the results are analyzed.

In the results portion of the paper we concentrate on three particular properties of bispectrum estimation. We show how the bispectrum estimate can be used to detect non-Gaussianity, non-linearity, and harmonic coupling. The detection of harmonic coupling is shown to be a useful property of the particular bispectrum estimation algorithm employed here.

14. Subject Terms.

signal processing, bispectrum and bicoherence of underwater acoustic signals, bispectral estimation algorithm

15. Number of Pages. 45

16. Price Code.

19. Security Classification of Abstract..

Unclassified

20. Limitation of Abstract.

NSN 7540-01-280-5500

Standard Form 298 (Rev. 2-89) Prescribed by ANSI Std. 239-18 298-102

Applications of Bispectral Analysis

A.M. Richardson W.S. Hodgkiss

Marine Physical Laboratory Scripps Institution of Oceanography San Diego, California 92152

Abstract

Estimates of the power spectral density have been found to be very useful in a variety of signal processing applications over the last several decades [1, 2, 3]. Higher order spectra contain information not present in the power spectrum and recently, estimates of higher order spectra have been shown to be useful in certain signal processing problems [4, 5, 6, 7, 8, 9, 10]. In particular, estimates of the bispectrum and bicoherence have been found useful in detecting non-Gaussianity and non-linearity [4, 6, 9], in system identification [11], and in detecting transient signals [10, 12].

This paper is concerned with the estimation of the bispectrum and bicoherence of underwater acoustic signals. We first introduce and describe some of the properties of the bispectrum and bico¬ herence. The bispectral estimation algorithm as implemented here at the Marine Physical Laboratory is then described in detail. Finally, several bispectra are estimated from actual data and the results are analyzed.

In the results portion of the paper we concentrate on three particular properties of bispectrum estimation. We show how the bispectrum estimate can be used to detect non-Gaussianity, non-linearity, and harmonic coupling. The detection of harmonic coupling is shown to be a useful property of the particular bispectrum estimation algorithm employed here.

Accesion For

NTIS CRA&I DT1C TAB Unannounced □

Justification _ _ _

Rv
Distribution f
Availability Codes
■	Avail and/or
Dist /H	Spe	cial

Contents

1 Introduction 1

1.1 Definitions . 1

1.1.1 Moments . 1

1.1.2 The Power Density Spectrum and the Bispectrum . 2

2 Estimation of the Bispectrum 3

2.1 Notation . 3

2.2 Welch’s Method of Power Spectrum Estimation . 3

2.3 Bispectrum Estimation . 4

2.4 The Bicoherence . 7

3 Properties of the Bispectrum and Bicoherence 8

3.1 Input/Output Relationships . 8

3.2 Symmetries in the Bispectrum Plane . . 8

3.3 Bispectra of Gaussian Processes . 8

3.4 Bicoherence of Linear Processes . 10

3.5 Bispectra of Harmonic and Quadratically Coupled Processes . 10

3.6 Statistical Properties of Bispectra . 11

4 Bispectral Analysis Results 18

4.1 Detection of a Non-Gaussian, Nonlinear Process . 18

4.2 Detection of a Harmonic Process . 18

5 Conclusions 42

6 Acknowledgements 43

*

1

1 Introduction

Power spectrum estimation has been a useful signal processing tool for many years [1, 2, 3]. Since the power spectrum is related to the autocorrelation it does not contain any information about moments higher than the second order. It is therefore of limited utility when one is trying to characterize or detect non- Gaussian noise. In system identification problems power spectral estimation techniques can only determine the magnitude response of the system. The phase response is unknown. (It is often assumed however that system is minimum-phase. If this assumption is correct then the complete response can be determined using power spectrum estimation methods.)

The bispectrum1 is a third-order spectrum which has a number of properties that potentially make it a valuable signal processing tool [13, 14, 15]. In practice the bispectrum has proven to be particularly useful in detecting non-Gaussianity and nonlinearity of ambient noise [4, 6, 7]. More recent work [10, 12] indicates that the bispectrum can also be useful in transient detection. The bispectrum can also be used to determine the phase response in system identification problems [11],

The remainder of this section is used to introduce notation and to define the bispectrum and bicoher¬ ence. In Section 2 we describe in detail the particular bispectral estimation algorithm implemented here at the Marine Physical Laboratory. Certain properties of the bispectrum and bicoherence are described in Section 3. We discuss a particular property of our bispectrum estimator which is useful in determining if spectral lines are harmonically related. We also discuss some of the statistical properties of the bispectrum. These properties are extremely useful in interpreting bispectral estimation results. In Section 4 we present results of bispectral estimation analysis as applied to data collected from a Swallow float and also to data collected from single elements of the MDA array.

1.1 Definitions

Let x(t) represent a real zero-mean stationary random process.

1.1.1 Moments

The first-order moment (or mean) is defined as

mx - £'{x(i)}, (1)

which, for the stated assumptions, is constant and equal to zero. The second-order moment (or autocor¬ relation) is

r**(r) = E{x(t)x(t + r)}. (2)

The bispectrum is related to the third-order moment (or autobicorrelation) which is defined as

Rxxx(n, r2) = E{x(t)x(t + Ti)x(t + t2)}. (3)

For the discrete sequence x[n] obtained by sampling the continuous-time signal x(t) every T seconds, i.e.

a:[n] = x(nT), the corresponding definitions are
mx =	£{*[«]}	(4)
r*®[n] =	£{a;[m]zr[m + n]}	(5)
RXxx[m, n] =	£{a:[/]x[/ + m]x[l + n]}.	(6)
'The function commonly called the bispectrum is more properly termed the bispectral density.

1

1.1.2 The Power Density Spectrum and the Bispectrum

The power density spectrum is defined as the Fourier transform of the autocorrelation function i.e.,

/oo

rxx(r)e~j2wfT dr. (7)

■CO

Assuming that x(t) is a voltage then the units of Pxx{f) are Volts2/Hertz (or assuming the power is dissipated in a 1 Ohm resistance, Watts/Hertz).

The bispectrum is defined as the Fourier transform of the autobicorrelation sequence,

/OO [OO

/ Rxxx(r1,T2)e-^^+^drldr2 (8)

■CO J — OO

with corresponding units of Volts3/Hertz2.

For discrete sequences the corresponding definitions are

PM = T £; rxx[n]e-^nT (9)

n=-oo

CO 00

JW/1./2) = E E fl***[™>«K;'2’(/,m+/3n)T- (10)

The value of the bispectrum depends on the amplitude of the component signals, and it is desirable to have a normalized form of the bispectrum [8]. The bicoherence is defined as [8]

bXXX (fuh)

_ Bxxxjfl 1 f 2) _

\/PM)PM)Pz*(fl+f2)‘

(11)

The bicoherence has units of Hertz-1/2

*

2

2 Estimation of the Bispectrum

Bispectral (and power spectral) estimation techniques can be classified as being either parametric or non- parametric (conventional) and also as being either direct or indirect [14]. Parametric estimators use a model (usually of the AR, MA, or ARMA type) for the process under investigation and obtain an estimate of the bispectrum by first estimating the parameters of the model. Non-parametric estimators make no assumptions about a model for the process. Parametric techniques can provide better estimates of the bispectrum than non-parametric techniques in situations where an assumption about the underlying model is valid [15] In particular, parametric estimation techniques are often used in the system identification problem. Direct type algorithms estimate the bispectrum directly from the data, indirect types compute the bispectrum from an estimate of the autobicorrelation.

The particular algorithm presented here falls into the non-parametric, direct estimation class of bis¬ pectrum estimators. It has many similarities to Welch’s method of power spectrum estimation, and so a brief overview of that technique is presented after a discussion of notation.

2.1 Notation

Let x[n\ denote the nth sample of the random process x(t), i.e.

a?[n] = x(nT) 0 < n < N — 1, (12)

where T is the sampling interval and it is assumed that we have N data points. In order to reduce the variance of a spectral estimate the data is usually divided into segments (unfortunately this also reduces the resolution). Assume that the data is divided into P segments of D samples with a shift of S samples between segments. The weighted pth segment can be written as

x(p)[n] = tu[n]x[n + p5] 0<n<D-l. (13)

The discrete-time Fourier transform (DTFT) of the weighted pth segment of data is equal to

o-i

x«(,)=r£ x^[n)e~j2vfnT. (14)

n=0

The factor of T is included in the definition so that the discrete-time Fourier transform will approximately equal the corresponding continuous-time transform. The discrete Fourier transform (DFT) of the same segment is equal to

o-i

X<*>[*] = Y, x^[n}e-^Tnk'D. (15)

n=G

The factor of T is omitted in the definition of the DFT to conform to standard notation.

2.2 Welch’s Method of Power Spectrum Estimation

There are several methods for estimating the power density spectrum. The various methods are discussed in detail in several texts on signal processing [1, 2, 3]. The bispectrum estimation procedure presented in the next section is analogous to Welch’s method for power density spectrum estimation and so the equations which define that method will be presented here. The notation closely follows that of Marple [1].

The sample spectrum of the weighted pth segment is given by

HPJ(f) = ^l*(p)(/)|2, (16)

3

where

(17)

D-l

Up — T ^2, U,2W-

n=0

The Welch estimate is obtained by averaging the sample spectrum from each segment

p-i

PrM) =4E

P—0

If a DFT is used to estimate the spectrum, the equations become

HPM = §V(P)MI2

p-i

p*m =

(18)

(19)

(20)

P—0

where X^p^[ifc] is the DFT of the pth weighted data segment.

For sinusoidal signals it is often more convenient to work with a function that indicates the power in each DFT bin. Denoting this function by Vxx[k] we have

vim =

p

p~ i

vm =

(21)

(22)

p= o

Vxx[k] has units of Volts2 (or Watts if measured across a 1 Ohm resistor).

2.3 Bispectrum Estimation

Our bispectrum estimation algorithm is based upon the following equation

1

Bxxx{fi,fi) = lim E{

where

M—*oo 1 (2 M + 1 )T M

XM(f) = T £ x[m]e-^nT.

(23)

(24)

m——M

The proof of this equation is fairly straightforward. After using the definition of X\i(f) in Eq.23 we have

M M M

Bxxx{fuh) = Jim

T2

M—co (2 M + 1) _

Y Y Y Rxxx[m-l,n-l\e-^^m-l^Te-^^n-^T. (25)

m=— M n=— JV/ /=— M

It can be shown that

M M M 2 M 2 M

X Y XI /[m -/,«-/]= X X (2Af+ 1-S[m,n])f[m,n], (26)

mzz-M n = -M l=-M

m=— 2M n — — 2M

where

S[m, n] = <

2 M + 1	m > 0; —2 M < n < —2 M + m
m — n	m > 0; —2 M + m < n < 0
m	m > 0; 0 <n <m
n	m > 0; m < n < 2 M.

(27)

4

The values of S[m, n] for m < 0 can be obtained from S[— m, — n] = S[m, n]. Note that S[m , n] is positive and is less than or equal to 2 M + 1. Using this in Eq. 25 gives

Bxxxifl 1/2) =

where it has been assumed that

CO 00

T. T S[m, n]Rxxx[m, n ] < 00. (31)

m=— 00 n = — 00

The last equation for the bispectrum is identical to Eq. 10 and therefore proves Eq. 23.

This proof suggests that the following algorithm for estimating the bispectrum. Let

1 P_1

JW/i - h) = p £ #&(/i - /a) (32)

p=0

where BxPxx(fi , h) is the sample bispectrum of the weighted pth data segment

B'ilifuh) = ~X^\h)X(p\f2)X^{h + f2) (33)

and X(p^(f ) is the discrete-time Fourier transform as defined previously in Eq. 14. In Eq. 23 t4 is equal to 1/(2 M + 1)T. It is difficult to say what t/j should be now because we are using windowed data that is finite in extent and our bispectrum estimate may be biased. The remainder of this section is concerned with finding the value of t/j which gives an unbiased estimator.

Assuming that E{BipJx{fi, /2)} does not depend on p (this is verified later) then the mean value of

BXxx(fi,h) is equal to

E{Bxxx(fuh)} = E{BipJx(fuf2)}. (34)

Substituting Eq. 33 yields

E{BXxx{h,h)} = ~E{X^{h)X^p\h)XW9{h + /2)} (35)

or (after using Eq. 14)

E{Bxxx{fi,h)} = 7T E E E E{x(P)[/]x(P)Hx^)[n]}e^'2,r^(m-0+A(n-01T- (36)

^ /=0 m= 0 n=0

The expectation on the right side of this equation is equal to

£{x^[/]x^[m]x^[n]} = RXXx{™ — l,n — /]u>[/]u;[m]u;[n]. (37)

Substituting into Eq. 36 gives

E{Bxxx{h, h)} = Tf (38)

6 ;=0 m=0 n=0

T 2

ZIYI ZIYI

lim , . . .

M-00 (2 M + 1) £

y y (2M + 1 - S[m, n])J?c„[m, n]e~^^m+^T (28)

m=-2M n = — 2A/

2 M 2 M cr ■«

2 V' /1 o[m, n]

lim T2 y y (1-

A/—00 ^ v

2M + 1

T2 y y Rxxxlm^le-WBm+h^T

)Rxxx[m, n]e-i2*V'">+hn)T

m=-2Af n=— 2A/ 00 00

(29)

(30)

m=— co n= — 00

5

The discrete-time third-moment function is the inverse Fourier transform of the discrete-time bispectrum

pl/(2T) ,1/(2 T)

Rx

/ Bxxx(Xu dXx dA2.

■1/(2 T) J — 1/(2T)

(39)

Substituting into Eq. 38 and simplifying gives (the limits of integration are omitted, but understood to be from —1/(2 T) to 1/(2 T) in the remainder)

E{Bxxx(f1,f2)} = ^J j Bxxx(\1,X2)W(f1-\1)W(h-\2)W'(f1+f2-\l-\2)d\1d\2 (40)

0—1

W(f) = T^2 w[n]e~j2rfnT ■ (41)

where

n= 0

This can be written more clearly as

E{Bxxx(fi,h)} = j-J j Bx„(\u\2)*(fi -A1,/2-A2)dAi d\2

(42)

where

*(/i,/a) = W{h)W{h)W*{h + /*)• (43)

The last equation for E{Bxxx(fx, f2)} indicates that the estimate is proportional to the true bispectrum convolved with $(/i,/2), i.e.,

E{Bxxx(fuh)} = Jj-BxxxUuh) * *(/l,/2),

(44)

and so our estimate is a biased estimator of the true bispectrum.

It is not possible to find a value of Ub that gives an unbiased estimate of the bispectrum for all frequency pairs. However, it is possible to find a value which gives an unbiased estimate of the mean cubed process. The third-moment sequence corresponding to our bispectrum estimate is

E{Rxxx[m,n]} = -^-Rxxx[m,n]'H[m,n] Ub

(45)

where 4([m, n] is the inverse Fourier transform of 4>(/i , f2). By setting Ub equal to '£[0,0] an unbiased estimate of #rrX[0,0] (the mean of the cubed process) is obtained. Ub (^[0,0]) can be determined from

Ub = / j$(fufi)dfidf2

= J J W(fx)W(f2)W*(fx + h) dfx df2

= T3 £ £ £ / e-i2*h{rn-i)T dfx f e-i**h(n-l)T df2

I ^ - n - n j

1=0 m=Q n= 0 0-1

= t£u,3[(],

(46)

(47)

(48)

(49)

1=0

the last expression being the most convenient.

If a DFT (see Eq. 15) is used to estimate the bispectrum, the corresponding equations are

B{£x[hM = ^-X^[kx}X^{k2}X^[kx + k2]

Ub

P-1

Bxxx[kx,k2] =

(50)

(51)

p= o

6

The three previous equations provide the basis for BISPCT, a computer program which estimates the bispectrum of a process from its sampled data [16].

When working with sinusoidal signals it is often convenient to define a scaled version of the estimated bispectrum that is equal to (see Eq. 22)

Bxcx[ki, £2] = ^jyp^2 ^2] • (52)

From the preceding equations we have

1 P_1

Bxxx[kuh] =

(53)

(54)

The units of Bxxx[ki, Jfe2] are Volts3.

2.4 The Bicoherence

We can estimate the bicoherence using estimates of the bispectrum and the power density spectrum,

(55)

The magnitude of the bicoherence is a convenient test statistic for use in the detection of coupled processes.

bXxx(fi ) f*i) —

_ Bxxxjfl I ft) _

f P**(h)P*x(h)Pxx(fl + h)

7

3 Properties of the Bispectrum and Bicoherence

3.1 Input/Output Relationships

If x(t) is the input to a linear system with impulse response h(t) and y(t) is the corresponding output, then the third-order moment of y(t) is related to that of x(t) by [11]

Ryyy(Tl > ^2)

-IS!

Rxxx{r 1 + t3 - ti,Ti + t3 - t3)h(ti)h(t2)h(t3) dti dt2<it3.

If we let H(f) represent the transfer function, the corresponding bispectra are related by [11]

Byyy{h,h) = Bxxx{h , f3)H (/l)tf (/2)iT ’(/t + /2)- The relation between the bicoherence functions is

byyyifu /z) — bxxx(fh /z)

H{h) H(h) £Ui±M \H(fl)\\H(f2)\\H(f1 + f2)\-

(56)

(57)

(58)

Note that the output bicoherence is affected only by the phase of the transfer function. (The magnitude of the bicoherence of the output is equal to the magnitude of the bicoherence of the input.)

3.2 Symmetries in the Bispectrum Plane

The third-order moment sequence of a stationary process has the following properties [11]

Rxxx(t1i T2) — Rxxx{T2i ri) = Rxxx(.~T 2, ~ 7"z) = Rxxx(~ fl , ~Tl + T2)

= Rxxx(-n + r2, -n) = RXXx(r 1 - r2, — r2). (59)

If the function Rxxx(ti , r2) is known in any of the six regions of Fig. 1 this set of equalities can be used to determine it everywhere. From the above symmetry relations and the definition of the bispectrum (Eq. 8) it follows that

BXxx(fi, h) — Bxxx(f2,fi) — Bxxx( fi — /2>/i) — Bxxx(—fi — f2,fi)

— Bxxx(f 2, fi /2) = Bxxx(fi, — fi — /2).

From Eq. 8 it also follows that the bispectrum obeys the following conjugate symmetry relation

Bxxx(fl,f2) = B*xxx(-f\,-f2)-

(60)

(61)

These relations imply that if the bispectrum is known in any of the twelve regions shown in Fig. 3.2 it can be determined everywhere in the plane.

From Eq. 10 it can be seen that the bispectrum of a discrete sequence is doubly periodic with period 1/T:

Bxxx(fi, /2) = Bxxx(fi + m/T, /2 + n/T). (62)

3.3 Bispectra of Gaussian Processes

The autobicorrelation of a stationary, Gaussian process is equal to zero which implies that the correspond¬ ing bispectrum is also zero. This property allows the bispectrum to be used as a tool for the detection of non-Gaussian processes. Note that a zero bispectrum does not imply a Gaussian process, i.e. processes which are not Gaussian may have zero bispectra. For example, any process whose samples are statistically independent and have a symmetric density function will have a zero autobicorrelation sequence and hence a zero bispectrum.

Hinich [17] has developed a test, based on the magnitude of the bicoherence estimate, for detecting non-Gaussianity.

8

3.4 Bicoherence of Linear Processes

A linear process is a time series that can be expressed in the form [6]

(63)

00

x[n ] = ^ a(m)e(n — m)

m= — oo

where e(n) is a purely random series, i.e. all of the {c(l), e(2), . . .} are mutually independent. Theoretically, the bicoherence of a linear time series is constant [7, 6].

A non-zero bispectrum (or bicoherence) implies non-Gaussianity. A non-constant bicoherence implies nonlinearity. It should be noted that this implies that all nonlinear processes are non-Gaussian.

Hinich [17] has developed a test, based on the spectral flatness of the bicoherence estimate, for detecting non-linearity.

3.5 Bispectra of Harmonic and Quadratically Coupled Processes

Assume that x[n] consists of a sum of three sinusoidal components,

3

x[n] = am cos(2 nfmTn + (64)

m=l

where /3 = f\ -f /2. The three component sinusoids could also be harmonically related (in which case all three frequency components will be multiples of a fundamental frequency). One situation in which the three components satisfy /3 = /j + /2 and are not necessarily harmonically related occurs in what is known as quadratic coupling. Quadratic coupling can occur when two time series (e.g. sinusoids with frequencies f\ and /2) pass through a transfer function that is nonlinear. Assuming we can approximate the input /output relationship by the first two terms in a Taylor series expansion as y(x) = ax + bx~ then, if the input consists of two sinusoids, the output will have frequency components at the original two frequencies and also at the sum and difference frequencies. The sum component at frequency /3 satisfies the above relationship. In addition, the phases satisfy the relation <t>3 — <f>i + <j> 2- This is termed quadratic phase coupling because it occurs whenever the output and input are related through a quadratic nonlinearity.

The discrete-time Fourier transform of the uniformly weighted pth segment of data, evaluated at frequency fi is approximately equal to

X(p)(/i) = a1DTej('t,l+2*hpTS) (65)

where S is the number of samples that we shift between overlapping segments. The corresponding expres¬ sions for X^p\f3) and X00(/3) are

A<p)(/2) = a2DTeH,h+2*f*pTS'> (66)

X(p)(/3) = a3DTej<'<t‘3+2,rf3pTS\ (67)

The estimated power density spectrum evaluated at these same frequencies is equal to (see Eqs. 16 and 18),

PM) = ai(DT) (68)

Pxx(h) = 4 (DT) (69)

Pxx(h) = a2(DT) (70)

(71)

The sample bispectrum (see Eq. 33) evaluated at ( f\ ,/2) is equal to

BxJx(fuh) = aia2a3(jDT)2e-'(?>1+<,!’3_^+2,r(/l+/2-/3)pT5). (72)

10

(73)

For harmonic and quadratically coupled processes in which fa = fa + fa, this reduces to

B(xpJx(fa,fa) = aia,MDT)2e^+^-^\

Note that the sample bispectrum in this case is independent of the segment index p.

If the three frequency components are a subset of lines from a harmonic process, it is usually assumed that the three phases are independent and identically distributed random variables uniformly distributed over a 2jt interval. In this case the ensemble average of the sample bispectrum is zero and therefore the bispectrum is also theoretically zero [13]. Fourth-order moments are then required to determine if the components are harmonically related [18]. However, our bispectrum estimator is based upon a time average of the sample bispectra instead of an ensemble average. Since the sample bispectrum does not depend upon the segment (time) index p, the value of our bispectrum estimator is

Bxxx(fa,fa) = (74)

The paper by Huber [8] recognizes the difference between the value of the estimator and the value of the true bispectrum in the case of harmonically related lines. They use the non-zero value of the bispectrum estimator as an aid in determining whether lines in a power density spectrum are harmonically related. The estimated bicoherence in this case is equal to (see Eq. 11)

bxxx(fi,fa) = VDTei<-*'+,h-*>\ (75)

It should be noted that in this particular problem our algorithm is a biased estimator of the true bispectrum at the point (fa, fa). Fortunately, this bias can be used to detect the presence of harmonics.

In the quadratic phase coupling problem the phases (j>\ and fa are statistically independent and fa = <j> i + fa. The phase of the sample bispectrum is equal to zero in this case and the ensemble and time averages of the sample bispectra are equal. The value of the bispectrum estimator in this case is

BXXX (fa, fa) = axdva-ziDT)2 , (76)

and the bicoherence estimate is equal to

bxxx(fa,fa) = y/DT. (77)

Points at which the magnitude of the bicoherence is equal to \fUf indicate the presence of either a harmonic or a quadratically coupled process.

This derivation assumes that the component frequencies are stable over the entire observation interval. If this is not the case, the bispectral and bicoherence peaks will be reduced in amplitude. Noise will also reduce the magnitude of the bicoherence.

3.6 Statistical Properties of Bispectra

As we have seen, a non-zero value in the bispectrum frequency plane is an indication of non-Gaussianity or frequency coupling. Of course, the presence of noise will also create peaks in the bispectral plane and so we must question whether or not a peak is statistically significant, i.e. indicative of non-Gaussianity or frequency coupling.

Hinich [17] provides a fairly complete statistical description of the bispectrum and bicoherence and describes in detail tests for non-Gaussianity and nonlinearity.

Huber, et al. [8], provide a simple approximation that is convenient for first-cut data analysis. Since bispectral estimates are asymptotically unbiased and asymptotically normal under mild conditions [19], they indicate that the variance of the off-diagonal elements of a bispectrum estimator of the type discussed here is equal to

D2T

var{B*xr(/i, fa)} = -jj-Pxx(fa)Pxx(fa)Pxx(fa + fa), (78)

11

Figure 4: Variance of the bispectrum of white Gaussian noise as a function of segment length (rectangular window). The top line and corresponding points (diamonds) are the theoretical and simulation results when no overlapping was used. The middle line and middle set of points (pluses) correspond to 50% overlap and the bottom line and points (squares) correspond to 75% overlap.

where Pxx(f) is the power density spectrum of the corresponding process. (The variance of a complex quantity 2 is defined here as var{z} = E{\z |2} — |i?{z}|2.) This result is applicable only for the case in which non-overlapping segments and rectangular windows are used. The variance of the diagonal elements is twice as large as that of the off-diagonal elements.

A computer simulation was performed to test the accuracy of the preceding theoretical result. The bispectrum of a white, Gaussian noise process of variance 0.25 Volts" was calculated. The sampling frequency is 1/2 Hertz. The corresponding power density spectrum is theoretically flat and equal to 0.5 Watts/Hertz. Rectangular windows and non-overlapping segments were used. According to the above equation the variance of the off-diagonal elements of the bispectrum is equal to

var {B,„(/i,/a)} = 2^. (79)

Figure 4 shows a graph of this equation along with the simulation results (shown as points) as a function of D for a value of N equal to 4096. Figure 5 graphs the variance as a function of N for a value of D equal to 256. (Only values of D and N that are equal to a power of two were used in the simulation.) There is good agreement between the approximate theoretical equation and the simulation results. The variance is also shown in the two figures for the cases in which the data windows were overlapped by 50% and 75%. To obtain the theoretical curves for these two cases the equation for variance with no overlap was multiplied by (1 — a/2) where a is the fractional overlap. (This is just a rule-of-thumb that seems to work well.) As expected the variance decreases as the amount of overlap is increased.

12

0)

o

G

(0

•H

<0

>

Figure 5: Variance of the bispectrum of white Gaussian noise as a function of observation length (rectan¬ gular window). The top, middle and bottom lines correspond to 0%, 50% and 75% overlap.

Huber also gives an approximate formula for the variance of the bicoherence [8]

var{6xxx(/i , /2)} = (80)

White Gaussian noise was used to test this equation for the same set of parameters used in the bispectrum simulation. The results are shown in Figures 6 and 7. When overlapping windows were used the same rule-of-thumb used to adjust the expression for the variance of the bispectrum was applied. Again there is very good agreement between the formula and the simulation results. Simulation also verified that the variance of the bicoherence estimate was relatively independent of the level of noise variance.

The data can be tapered at the ends of the window to further reduce the variance. Figs. 8 through 11 show simulation results when a Kaiser-Bessel window with parameter /? = 2.5t is used. These figures should be compared to Figs. 4 through 7.

Since the bicoherence is asymptotically normal, |6xxx(/i , /2)|2 will be asymptotically \2 distributed with two degrees of freedom for a process with a vanishing true bispectrum (for example, a Gaussian process). Since the expected value of \bxxx(f\ , /a))2 is equal to var{&xxx(/i , /2)} in this case the distribution is known and we can determine levels at which the magnitude squared bicoherence is statistically significant. For example, the probability that \bxxx(f\, /2)|2 exceeds three times its mean value is less than 0.05, i.e.

Prob{j6xxx(/i, /2)|2 > 3var{6xxx(/i, /2)}} < 0.05 (81)

or

Prob{|6xxx(/i,/2)| > >/3var{6xxx(/i,/2)}} < 0.05. (82)

In order to detect quadratic or harmonic coupling it is necessary for the significance level to be less than sfUT (see Eqs. 75 and 77). As a rule-of-thumb we select the bispectrum estimation algorithm parameters

13

2500

2000

u 1500

c

<U

u

g 1000

500

0

0 500 1000 1500 2000

D

Figure 6: Variance of the bicoherence of white Gaussian noise as a function of segment length (rectangular window). The top, middle and bottom lines correspond to 0%, 50% and 75% overlap.

<D

U

C

<0

•H

S-l

m

>

800	- 1 - 1 - 1 -	- , -
700	i	■
600	•1	-
500	■9 ti]	-
400		-
300	\\ \\	-
200	T*	■
100 n	A\ \v\	-

0 4000 8000 12000 16000

N

Figure 7: Variance of the bicoherence of white Gaussian noise as a function of observation length (rectan¬ gular window). The top, middle and bottom lines correspond to 0%, 50% and 75% overlap.

14

Figure 8: Variance of the bispectrum of white Gaussian noise as a function of segment length (Kaiser- Bessel window). The top set of points (diamonds) are simulation results when no overlapping was used. The middle set of points (pluses) correspond to 50% overlap and the bottom points (squares) correspond to 75% overlap.

15

40

<1>

o

G

<0

•H

U

>

35

30

25

20

15

10

5

0

□	- 1 -	- , -	1
“+				-
o
0				-
- +				•
o o
a •				-
£	o £		» -
0	4000	8000	12000	16000

N

Figure 9: Variance of the bispectrum of white Gaussian noise as a function of observation length (Kaiser- Bessel window). The top, middle and bottom lines correspond to 0%, 50% and 75% overlap.

O

G

m

•H

U

<0

>

1000

900

800

700

600

500

400

300

200

100

0

-		i	i 1	- 1 — 0
-				+
-				□ ’
-			o	-
			+
-			□	-
		o
	ffl	_ g _	1 . 1	■ . , n- I. ..i.n

0 500 1000 1500 2000

D

Figure 10: Variance of the bicoherence of white Gaussian noise as a function of segment length (Kaiser- Bessel window). The top, middle and bottom lines correspond to 0%, 50% and 75% overlap.

16

N

Figure 11: Variance of the bicoherence of white Gaussian noise as a function of observation length (Kaiser- Bessel window). The top, middle and bottom lines correspond to 0%, 50% and 75% overlap.

( N , D, etc.) such that the variance in the white Gaussian noise case satisfies

v/3var{6 XXX (/i,/a)} < v'DT/ 5

var{6xxx(/i,/2)} < DT/75. (84)

From Eq. 80 this implies, that for non-overlapping, rectangular data windows D and N should satisfy

# * W <85>

For practical reasons it may be necessary to use overlapping and nonrectangular windows in order to increase the ratio of D to N. In that case, simulation is used to determine the estimation algorithm parameters such that Eq. 84 is satisfied.

17

4 Bispectral Analysis Results

We now present results from applying bispectral analysis to real ocean acoustic data. In each case, graphs of the power density spectrum, the magnitude of the bispectrum, the bispectrum normalization function, and the magnitude of the bicoherence are presented. The bispectrum normalization function is the denominator of Eq. 11. Graphs of this function are useful in interpreting bicoherence results. The power density spectrum is plotted in decibels. The bispectrum and the bispectrum normalization are both plotted on the same scale. We plot 20/3 times the base 10 logarithm of both functions and the peak of the bispectrum normalization function is defined as 0 dB. (The 20/3 value is just a convenient factor to use for comparing bispectrum and power density spectrum levels.) The magnitude of the bicoherence is normalized by dividing by \/DT the result is then plotted on a linear scale with range between 0 and 1.

4.1 Detection of a Non-Gaussian, Nonlinear Process

The bispectrum of a data set from the August 1988 Swallow float experiment [20] was calculated. These data were recorded by float 5 which was at a depth of about 700 m. The bottom depth was approximately 3800 m. The Swallow floats contain three geophones and one hydrophone. Only the bispectrum of the hydrophone is presented here. The Swallow floats use a sampling rate of 50 Hz. A single record (2250 samples or 45 seconds of data) was selected for bispectral analysis. To satisfy Eq. 84 in the white Gaussian noise case, Eq. 85 suggests a segment length of only 30 points if non-overlapping, rectangular windows are used. Simulation indicates that the variance of the bicoherence magnitudes of white Gaussian noise with N = 2240, D = 128, 50% overlap and a Kaiser-Bessel window is equal to 0.035. This satisfies Eq. 84 and substituting into Eq. 82 gives

Prob{|ixxx(./i,/2)|> 0.32} < 0.05, (86)

or, for the normalized bicoherence value,

Prob{\bxxx(flth)\/VDT > 0.2} < 0.05 (87)

The results are shown in Figs. 12 through 19. Surface and gray-scale plots are shown for the estimated bispectrum, bispectrum normalization and normalized bicoherence. Two surface plots of the normalized bicoherence are shown. The second surface plot has a raised “floor” at 0.2. Points in the normalized bicoherence above the floor are statistically significant at the 0.05 level. In the normalized bicoherence gray-scale plots, values at or above the 0.2 level are shown as white.

There are approximately 1024 points in the principle domain of the bicoherence. If the data were white Gaussian2 we would expect approximately 50 points in our bicoherence estimate to have values above 0.2 (which corresponds to the 0.05 statistical significance level). There are actually 124 points that have values greater than 0.2. We only expect five points to exceed 0.27 (which corresponds to the 0.005 significance level). There are 58 points in our estimated bicoherence that exceed this level. Examination of the bicoherence estimates shown in Figs. 17 through 19 indicates that the low-frequency portion of the signal (below 5 Hz) is non-Gaussian and nonlinear. The Hinich tests [6] could be applied to further aid in deciding whether or not the time series is in fact non-Gaussian and nonlinear.

For comparison the power density spectrum, bispectrum, and bicoherence of the next data record was also calculated. The results are shown in Figs. 20 through 27. The bicoherence estimate for this data record looks much like the previous one. There are more points that exceed the 0.2 value (152), but a fewer that exceed the 0.27 value (54).

4.2 Detection of a Harmonic Process

A power density spectrum for a particular set of data taken from the 1991 MDA experiment is shown in Fig. 28. (This particular data set was obtained from element 17 of leg A of the array on day 196 starting at

2 If the data were correlated Gaussian we would expect even fewer values to exceed the 0.2 level.

18

Figure 12: Power density spectrum of the acoustic pressure as measured by a Swallow float. First data record.

19

fx (Hz)

Magnitude (dB)

(Hz)

23

Figure 20: Power density spectrum of the acoustic pressure as measured by a Swallow float. Second data record.

24

Frequency (Hz)

Figure 28: Power density spectrum of the acoustic pressure as measured by element 17 of leg A of the MDA array.

time 00:52:39.483.) This data set was selected because of the set of lines that appear in the power density spectrum at multiples of approximately 9.1 Hz. Their appearance was correlated with the sighting of a large ship transiting the area. The power density spectrum was estimated using Welch’s method and a 512 point FFT. The data was tapered (using a Kaiser-Bessel window with parameter f3 — 2.5t) and averaged using 50% overlap. The sampling frequency is 150 Hz. Using these same parameters, simulation indicates that the variance of the bicoherence of white Gaussian noise is equal to 0.044 for an observation time of 64 seconds (9600 samples). This satisfies Eq. 84 and substituting into Eq. 82 gives

Prob {\bxtx(fi,f2)\ > 0.363} < 0.05 (88)

or

Prob{|ixxx(/1,/2)|/v/£T> 0.2} < 0.05. (89)

The bispectrum and the bispectrum normalization function (the denominator of the bicoherence) for this data set are shown in Figs. 29 through 32. Examination of the power density spectrum indicates that a harmonic process with a fundamental frequency of 9.1 Hz may be present. There is a plainly visible grid in Fig. 30. The regularity of the grid spacing provides an indication that the process contains many harmonically related components.

The normalized bicoherence for this case is shown in Figs. 33 through 35. In Fig. 34 there are several peaks above the 0.2 bicoherence level. This is to be expected because of the large number of points in the principle domain of the bispectrum plane. For an FFT length of 512 points there are 16512 (N(N +4)/16) points in the principle domain of the bispectrum plane. For a white Gaussian noise process we expect

29

(Hz)

70 0

Figure 29: Bispectrum of the acoustic pressure as measured by element 17 of leg A of the MDA array.

f* (Hz)

Figure 30: Bispectrum of the acoustic pressure as measured by element 17 of leg A of the MDA array.

30

Magnitude (dB)

Figure 34: Bicoherence of the acoustic pressure as measured by element 17 of leg A of the MDA array. Floor raised to 0.2.

0.10

Figure 35: Bicoherence of the acoustic pressure as measured by element 17 of leg A of the MDA array

Peak	Value	h	h
1	0.4671600	9.0820312	9.0820312
2	0.4640700	36.6210938	9.3750000
3	0.4503400	18.4570312	9.0820312
4	0.4412500	27.5390625	27.2460938
5	0.4384600	34.8632812	9.0820312
6	0.4345900	39.8437500	9.0820312
7	0.4082000	29.8828125	9.0820312
8	0.4057400	37.5000000	9.3750000
9	0.3947200	54.4921875	9.0820312
10	0.3864100	18.1640625	14.6484375
11	0.3838200	27.2460938	18.1640625
12	0.3713800	23.7304688	18.1640625
13	0.3636300	40.1367188	28.1250000
14	0.3618500	38.3789062	9.0820312
15	0.3472300	31.0546875	29.5898438
16	0.3445000	40.4296875	0.2929688
17	0.3422700	44.2382812	9.0820312
18	0.3416200	45.9960938	18.1640625
19	0.3414100	33.9843750	13.7695312
20	0.3396300	54.7851562	11.1328125

Table 1: The 20 largest peaks in the bicoherence of the acoustic pressure as measured by element 17 of leg A of the MDA array.

826 points to exceed the 0.2 bicoherence value. There are actually 1273 points in Fig. 34 that exceed the 0.2 value. This suggests that the data may have some non-Gaussian characteristics. The 20 largest peaks in the bicoherence and their locations are listed in Table 1. The four largest peaks all occur within a single FFT bin width of multiples of 9.082 Hz. The largest peak has coordinates (9.082, 9.082) in the /i, fi frequency plane. This is an indication of coupling between the 9.082 and 18.164 Hz components. The second largest peak in the bicoherence occurs at (36.621, 9.375). The 9.375 value is one FFT bin away from the 9.082 value and the 36.621 value is one FFT bin away from four times 9.082. This indicates coupling between the fundamental, the third, and fourth harmonic components. The third highest peak occurs at (18.457, 9.082) and indicates coupling between the fundamental, the first, and the second harmonics. The fourth peak is located at (27.539, 27.246) and implies coupling between the second and fifth harmonics. The next two peaks occur in the area about (36.621, 9.082) which is in the region of the third harmonic. In Fig. 34 we can see considerable spread about the third harmonic along the /a = 9.082 line. Referring to the plot of the power density spectrum we can see that the third harmonic is much weaker than the others. Going further down the table one can find other peaks that are located near multiples of 9.082 Hz. The ninth peak, for example, occurs at (54.492, 9.082). The value 54.492 is exactly six times greater than 9.082 and indicates coupling between the fundamental and the fifth and sixth harmonics. Most of the time, if a peak does not have coordinates which are exact multiples of 9.082, it is located in a bin with a slightly higher frequency. This indicates that the fundamental frequency is slightly higher than 9.082 Hz.

As another example, data collected at the same time but from a different element of the array (element 33) were analyzed using the same parameters as in the previous case. The results are shown in Figs. 36 through 43.

There are 1495 points that exceed the 0.2 bicoherence value for this element. Again, if the process were white and Gaussian we only expect 826 points to have values above 0.2. Looking at the table of

34

Figure 36: Power density spectrum of the acoustic pressure as measured by element 33 of leg A of the MDA array.

35

Magnitude

Magnitude (dB)

Figure 42: Bicoherence of the acoustic pressure as measured by element 33 of leg A of the MDA array Floor raised to 0.2.

0.15

0.10

Figure 43: Bicoherence of the acoustic pressure as measured by element 33 of leg A of the MDA array.

Ltude

Peak	Value	fi	h
1	0.5725200	9.3750000	9.0820312
2	0.4942200	18.4570312	9.0820312
3	0.4583200	19.3359375	12.8906250
4	0.4269600	31.6406250	6.1523438
5	0.4207500	54.7851562	9.0820312
6	0.4202700	0.2929688	0.2929688
7	0.3956600	27.2460938	9.3750000
8	0.3853800	45.4101562	9.0820312
9	0.3849600	19.0429688	13.4765625
10	0.3826200	20.8007812	11.4257812
11	0.3709000	48.6328125	0.2929688
12	0.3685600	13.1835938	12.8906250
13	0.3659900	34.5703125	17.8710938
14	0.3571500	36.6210938	36.3281250
15	0.3545500	53.3203125	15.8203125
16	0.3530600	18.1640625	9.6679688
17	0.3527000	13.1835938	10.5468750
18	0.3519800	57.1289062	9.0820312
19	0.3518800	37.5000000	31.6406250
20	0.3511600	16.1132812	15.5273438

Table 2: The 20 largest peaks in the bicoherence of the acoustic pressure as measured by element 33 of leg A of the MDA array.

40

peak values and their locations we can again see that several of the peaks have coordinates that are near multiples of 9.082. This provides further evidence that component lines in the power density spectrum are harmonically related.

41

5 Conclusions

Bispectral analysis appears to be a promising technique for use in detecting non-Gaussianity and non¬ linearity. Although the bispectrum of a harmonic process is zero the particular bispectrum estimation algorithm discussed here has certain properties that make it useful in detecting harmonics.

The bicoherence plays a key role in the detection of non-Gaussian and non-linear characteristics and also in detecting frequency coupling. Knowledge of the variance of the estimated bicoherence under the Gaussian hypothesis is required to determine if peaks in the bicoherence are statistically significant. For the estimation algorithm discussed here, a simple expression for the variance exists when rectangular windows are used. For non-rectangular windows, a simulation is required to determine the variance.

42

6 Acknowledgements

The authors are indebted to Gerald D’Spain and Jui-Yuan Chang for several useful discussions and support during the course of this work. This work was supported by the Office of Naval Technology under NRL contract #N000 14-92- K-2005.

43

References

[1] S. Lawrence Marple, Jr. Digital Spectral Analysis with Applications. Prentice-Hall, Englewood Cliffs, NJ, 1987.

[2] Steven M. Kay. Modem Spectral Estimation. Prentice-Hall, Englewood Cliffs, NJ, 1988.

[3] John G. Proakis, Charles M. Rader, Fuyun Ling, and Chrysostomos L. Nikias. Advanced Digital Signal Processing. Macmillan, New York, 1992.

[4] Klaus Hasselmann, Walter Munk, and Gordon MacDonald. Bispectra of ocean waves. In Murray Rosenblatt, editor, Proceedings of the Symposium on Time Series Analysis, pages 125-139. John Wiley and Sons, Inc., June 1962.

[5] M.J. Hinich and C.S. Clay. The application of the discrete fourier transform in the estimation of power spectra, coherence, and bispectra of geophysical data. Reviews of Geophysics , 6(3):347-363, August 1968.

[6] Patrick L. Brockett, Melvin Hinich, and Gary R. Wilson. Nonlinear and non-gaussian ocean noise. J. Acoust. Soc. Am., 82(4) : 1386— 1394, October 1987.

[7] Melvin J. Hinich, Davide Marandino, and Edmund J. Sullivan. Bispectrum of ship-radiated noise. J. Acoust. Soc. Am., 85(4):1512-1517, April 1989.

[8] Peter J. Huber, B. Kleiner, T. Gasser, and G. Dumermuth. Statistical methods for investigating phase relations in stationary stochastic processes. IEEE Transactions on Audio and Electroacoustics, AU-19(l):78-86, March 1971.

[9] Melvin J. Hinich and Gary R. Wilson. Detection of non-gaussian signal in non-gaussian noise using the bispectrum. IEEE Transactions on Acoustics, Speech and Signal Processing, 38(7):1126— 1131, July 1990.

[10] Melvin J. Hinich. Detecting a transient signal by bispectral analysis. IEEE Transactions on Acoustics, Speech and Signal Processing, 38(7):1277— 1283, July 1990.

[11] Athanasios Papoulis. Probability, Random Variables, and Stochastic Processes. McGraw-Hill, third edition, 1991.

[12] Lisa A. Pflug, George E. Ioup, Juliette W. Ioup, and Robert L. Field. Properties of higher-order correlations and spectra for bandlimited, deterministic transients. J. Acoust. Soc. Am., 91(2):975— 988, February 1992.

[13] Jerry M. Mendel. Tutorial of higher-order statistics (spectra) in signal processing and system theory: Theoretical results and some applications. Proceedings of the IEEE, 79(3):278-305, March 1991.

[14] Chrysostomos L. Nikias and Mysore R. Raghuveer. Bispectrum estimation: A digital signal processing framework. Proceedings of the IEEE, 75(7):869— 891 , July 1987.

[15] Mysore R. Raghuveer and Chrysostomos L. Nikias. Bispectrum estimation: A parametric approach. IEEE Transactions on Acoustics, Speech and Signal Processing, ASSP-33(4):1213-1230, October 1985.

[16] Jui-Yuan Chang. BISPCT. The computer reference manual page for the BISPCT program, Scripps Institution of Oceanography, Marine Physical Laboratory.

[17] M.J. Hinich. Testing for gaussianity and linearity of a stationary time series. J. Time Series Anal., 3(3): 145— 158, 1980.

44

[18] Chrysostomos L. Nikias and Jerry M. Mendel. Signal processing with higher-order spectra. IEEE Signal Processing Magazine, 10(3): 10— 37, July 1993.

[19] D.R. Brillinger and M. Rosenblatt. Asymptotic theory of estimates of 1-th order spectra. In B. Harris, editor, Advanced Seminar on Spectral Analysis of Time Series, pages 153-188. Wiley, New York, 1967.

[20] G.L. D’Spain, R.L. Culver, W.S. Hodgkiss, and G.L. Edwards. Freely drifting swallow float ar¬ ray: August 1988 trip report. Technical report, Marine Physical Laboratory, Scripps Institution of Oceanography, San Diego, CA, 1989. MPL Tech. Mem. 420.

[21] K.S. Lii and K.N. Helland. Cross-bispectrum computation and variance estimation. ACM Transac¬ tions of Mathematical Software, 7(3):284-294, September 1981.

45

ONR/MPL Report Distribution

Administrative Grants Officer (1)

Office of Nava! Research Resident Representative N6601 8 University of California, San Diego (Mail Code 0234) 8603 La Jolla Shores Drive San Diego, CA 92093-0234

Director (2)

Naval Research Laboratory Atten: Code 2627, 5160 Washington, D.C. 20375

Defense Technical Information Center (4) Building 5, Cameron Station Alexandria, VA 22314